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1 INTRODUCTION 



Q , ABSTRACT 

■ We present a theory for describing the evolution of a galaxy caused by stochastic 
| events such as weak mergers, transient spiral structure, orbiting blobs, etc. This noise 

■ excites large-scale patterns that drives the evolution of the galactic density profile. 

In dark-matter haloes, the repeated stochastic perturbations preferentially ring the 
lowest-order modes of the halo with only a very weak dependence on the details of 
their source. Shaped by these modes, the profile quickly takes on a nearly self-similar 
form. We show that this form has the features of the "universal profile" reported by 
Navarro, Frenk, & White independent of initial conditions in a companion paper. In 
this sense, this noise-driven process is a near-equilibrium form of violent relaxation. 

Key words: galaxies:evolution — galaxies: haloes — galaxies: kinematics and dy- 
namics — cosmology: theory — dark matter 
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' Both numerical experiments and observations suggest a universality in galaxy profiles resulting from dissipationless collapse. 
(-H | Remnants from merger simulations (e.g. van Albada 1982, Barnes 1989) typically manifest the often observed r x ^ 4 -type profile 
O i 1 in elliptical galaxies. More recently, the same feature has arisen in the context of cosmological large-scale structure simulations. 
Researchers beginning with Navarro, Frenk & White (1997) have noted that CDM collapse results in a haloes with a similar 
overall profile. Debate remains about the precise functional form of these profiles, but these and the r x ' 4 -law are all quite 
similar, suggesting that some general violent-relaxation-like dynamical mechanism is at work. 

A precise explanation of the dynamics of violent relaxation, the process driving the convergent evolution in by large 
fluctuations in potential, remains elusive. There have been several attacks on this problem. Tremaine, Henon & Lynden-Bell 
(1986) follow a statistical mechanics approach and explore extremizing functionals (see their paper for a review of prior work). 
Spergel & Hernquist (1992) extended this approach with the additional ansatz that the orbit-perturbation interaction can 
be approximated by kicks near perigalacticon. Recently, this has been followed by Mangalam, Nityananda & Sridhar (1999) 
' who incorporate these physical mechanisms in a semi-analytic evolution equation based on diffusion in action space. Over the 
same period, Kandrup has pursued a mathematically rigorous exploration into the geometry of phase space resulting from 
the details of Hamiltonian flows (e.g. see Kandrup 1998 and references therein). 

In this paper, I take a different approach to the problem. Rather than trying to identify an appropriate form for the 
phase-space distribution to start, I treat the problem as an initial value problem for a stochastic process and develop an 
evolution equation for exploring its consequence on stellar systems. A companion paper (Paper 2) will apply this to the 
evolution of galactic haloes in particular. The underlying motivation is as follows. Previous work by the author has shown 
that fluctuations in stellar systems on the largest scales can be strongly amplified by their own self gravity (e.g. Weinberg 
1993). This means that large-scale fluctuations will greatly exceed their Poisson amplitudes. In particular, Weinberg (1993) 
explored the idealized case of periodic cube; the fluctuations become very large as the system size approaches its Jeans' 
length. Similarly for stable galaxies, the fluctuations in a system will be largest at its discrete modes. In addition, Weinberg 
(1994) argues that galaxies will often have very weakly damped m = 1 (sloshing or seich) modes and these result in large 
excitations when excited (see Vesperini & Weinberg 2000 for another example). Putting this together, one might ask: if noise 
preferentially excites particular modes independent of the noise source, is it possible the repetitive stochastic response of the 
galaxy will lead to some characteristic features, independent of its initial conditions? 

In order to answer this question, this paper concentrates on the theoretical framework for describing the evolution due 
to stochastic fluctuations. Beginning with a description of the linear response of a galaxy to excitation, and assuming that 
the process is Markovian, one may expand the Boltzmann collision term in a series. Analogous to two-body collisions, only 
the first two terms contribute and the resulting evolution equation has a Fokker-Planck form. This equation is very far from 
being analytically tractable because the diffusion coefficients depend on integrals over all of phase space under the stochastic 
perturbation. It does not solve the classic violent relaxation problem because this approach is perturbative and assumes 
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that the stellar system remains near an equilibrium, Nonetheless it does incorporate the same underlying processes, the self- 
gravitating response to noise, and yields profiles with many of the same features that are found in the numerical experiments 
as Paper 2 will demonstrate. 

The plan for this paper is as follows. The overall approach is outlined in §^l] followed by an explicit derivation for 
spherical haloes in £2.2. This could be straightforwardly extended to many standard geometries geometry (e.g. disk or disk 
and halo together). We describe the character of the noise for two general cases, transient and quasi-periodic perturbers (a 
halo of black holes, for example) in These two cases result in qualitatively different behavior and represent the most 
plausible astronomical scenarios. Transient noise is probably the most relevant and important. Finally, the main features are 
summarized in Paper 2 will review the basic physics and apply these methods to understanding the evolution of halo 
in a noisy environment, e.g. just after formation, and may be a better place to begin for those interested in astronomical 
consequences rather than the kinetic theory. 



2 DERIVATION OF THE EVOLUTION EQUATION 
2.1 Overview 

The outer halo in large galaxies is less than 10 dynamical times old so primordial inhomogeneities will not have had time 
to phase mix and continued disturbances from mergers will not have relaxed (see Tremaine 1992 for additional discussion). 
These, as well as any intrinsic sources of noise, e.g. a population of 10 6 Mq massive black holes, dwarf galaxies, debris streams 
(Johnston 1998, Morrison et al. 2000) or dark clusters, are amplified by the self-gravity of the halo. In the current epoch, these 
distortions create potentially observable asymmetries in the stellar and gaseous Galactic disk. Just after galaxy formation, 
this noise may be sufficient to drive the evolution of the halo (Paper 2). Other possible applications include the evolution of 
proto-stellar clusters and the evolution of the proto-stellar binary distribution in the noisy star forming environment (work in 
progress). All of these problems motivated the development of a stellar dynamical kinetic equation that can handle general 
stochastic processes. 

The general problem will be familiar to most dynamicists and could be solved following the standard two-body approach: 
beginning with the collisional Boltzmann equation, one writes the collision term in Master equation form and expands in 
a Taylor series to derive Fokker-Planck equation (e.g. Binney & Tremaine 1987, Spitzer 1987). For a spherically symmetric 
system, the phase-space distribution is a function of two action and and two angle variables. Averaging over times short 
compared to the relaxation time but long compared to the dynamical time (orbit-averaging), one obtains a Fokker-Planck 
equation. 

Alternatively, recent work in statistical mechanics and noise theory in general has developed a body of methods for treating 
stochastic differential equations based directly on transition probabilities. If P(x',t + r\x,t) is the transition probability to 
some new state x at time t + r from the initial state x at time t, then the following integral equation determines all subsequent 
evolution of the distribution f(x): 

f(x,t + r) = / dx'P(x,t + t\x' ,t)f(x' ,t). 



By expanding the transition probability in its moments of x — x' for small r, one may derive a differential equation of the 
form: 

n = l 

where the coefficients 7?'™' are proportional to the time-derivative of the moments in transition probability. This is known as 
the Kramers-Moyal expansion. 

For galaxy evolution, the state variable x in the equations above is replaced with the phase-space six vector. Since we 
are interested in long-term evolution, we may orbit average the evolution equation. Writing the phase space as actions and 
angles turns the orbit average into an angle average. For any given noise process, we can solve for the change in actions of any 
orbit in the galaxy using a perturbative approach and similarly derive the change in the phase-space distribution function as 
described in Weinberg (1998). Evaluation of these quantities lead directly to the moments needed for the Fokker-Planck- type 
evolution equation. The most subtle aspect of the development below is enforcing a consistent time ordering. Implicit in the 
angle averaging is a short dynamical time scale and a long evolutionary time scale. The stochastic perturbations are on short 
time scales and instantaneous from the evolutionary point of view. Although this approximation may seem restrictive and 
only marginally true for some scenarios of interest, it yields results in good agreement with n-body simulation in the case of 
globular cluster evolution and the closely related case of self-gravitating fluctuations explored in Weinberg (1998). 

After deriving the evolution equation below, we work out the details of the Fokker-Planck coefficients for two examples 



m 



§^: transient distortions such as fly-by encounters and quasi-periodic distortions such as orbiting super-massive black holes. 



The consequences of both of these will be explored in Paper 2. 
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2.2 Derivation of the Fokker-Planck equation 

The main advantage of the Kramers-Moyal expansion is that the noise process appears explicitly in an initial-value form. 
Because the response of the galaxy to a stochastic event can be straightforwardly computed in perturbation theory or by 
n-body simulation, we can compute the noise-driven evolution for a wide variety of scenarios. For this reason, it is better 
suited to treating general stochastic noise than the Master equation, even though the two approaches are formally equivalent. 

Although the Kramers-Moyal expansion has an infinite number of terms in general (cf. eq. |TJ), the Pawula Theorem 
(Pawula 1967) shows that consistency demands that the expansion either stops after two terms and takes the standard 
Fokker-Planck form or must have an infinite number of terms. For the stellar dynamical case considered here, the series 
truncates after two term s. The approach approach sketched here is clearly described by Risken (1989, Chap. 4). 



As outlined in 



2.1 



one begins the derivation with the definition of transition probability: 

/(I, t + t) = Jdl' P(I, t + t\I', t)f(f,t) (2) 

where an average over the rapidly oscillating angles is implied and P is the conditional probability that a state has I at time 
t + t if it has I' at time t initially. A Taylor series expansion of the integrand in A = I' — I followed by a change of variables 
and integration over A. In the limit r — > this expansion leads directly to 



0/(1, t + r) v^/ d 



dt dl 

71 = 1 



(I, *)/(!,*). (3) 



Note that P is the probability of a change in I due to stochastic events. Therefore, the formal time derivatives in the expansion 
and phase space integral is better considered as the limit for small r (but for r greater than a dynamical time) of the ensemble 
average of stochastic events. If we let £ be the stochastic value of I, the expansion coefficients describing the stochastic variables 
€ is: 



D^(x,t)= i limi ([£(* + r)- I] 



In TV dimensions, this becomes: 



(4) 

«*)=! 



where the moments M'™' are 

M i"j a ....,j»( x »*» T ) = (( J ii - - ij*)- ■ ■ ih n - Ijj) ■ (6) 

The angle brackets denote the integration over the conditional probability of obtaining the state variable Ij at time t + r 
given Ij at time t. 

We will see below that the response function guarantees the that the Kramers-Moyal expansion terminates after two 
terms in the limit of weak perturbations and will assume this here. This is consistent with our intuition that repetitive weak, 
stochastic excitation in a galaxy is a Markov process. The evolution equation is then the Fokker-Planck equation in standard 
form: 

*^M=L, P (M)/(M) (7) 
where 

L„(I,*) = -A J? « (I>t) + _g_ I) g) (I>t) . (8) 

We see that the natural coordinates for the Boltzmann equation are action-angle variables. For a collisionless equilibrium, 
the actions are constant and the angles advance at constant rate. In deriving this Fokker-Planck equation, one averages over 
orbital periods in favor of following the evolution over longer time scales. The distribution function is then a function of 
actions alone, / = /(I). 

However, we will concentrate on the evolution of haloes modeled as a collisionless spherical distribution for application 
in Paper 2 and adopt the traditional E,J,J Z or E, k = J/Jmax(E), cos /3 = J z /J phase-space variables for computational 
purposes. In addition, we will not consider any processes with a preferred axis so we can average equation (Q) over (3 with no 
loss of information. If we are willing to restrict ourselves to an isotropic distribution, we can average over k to yield a 1 + 1 



dimensional Fokker-Planck equation in time and energy, E\ this is described in §2.3 below. Because the diffusion coefficients 
depend on the distribution function, the Fokker-Planck equation in ([?]) is non-linear, just as for globular cluster evolution. 
However, it is straightforward to solve this numerically by iteration. 

To derive the coefficients and D^ 2 \ we will use the method described in Weinberg (1998, hereafter Paper 1). To 

summarize, we represent distortions in the structure of halo in a biorthogonal basis. Any distortion can be summarized then 
by a set of coefficients in three indices. Because large spatial scales are most important in understanding global evolution, we 
can truncate this expansion and still recover most of the power. Moreover, we can analytically compute the self-gravitating 
response of the halo to some arbitrary perturbation as previously described. This development gives us f (t) for all phase-space 
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variables (cf. eq. ^) and the appropriate ensemble averages give us the required diffusion coefficients. For example, if one uses 
the same biorthogonal basis in an n-body simulation of a desired transient process, the time series of coefficients can be used 
directly to derive the coefficients D^ 1 ' and I)' 2 ' after removing the time-invariant (DC) component which corresponds to the 
equilibrium background, 
substructure, 



We will consider point mass perturbers in 
spreading debris trails) in 



3.2 



i.l and transient perturbers (dwarf galaxies, decaying 



2.3 Averaged Fokker-Planck equation 



A number of authors have described transformation of the multivariate Fokker-Planck equation (Rosenbluth et al. f957, 
Risken 1989). The approach is the familiar one: write the equation in terms of scalars, covariant and contravariant vectors 
and tensors and covariant derivatives only. In the first case, the authors use the Jacobian of the coordinate transformation as 
a metric and in the second, the authors use the diffusion matrix. We will use the first case here. Denote the Jacobian of the 
coordinate transformation as J. Under a change of coordinates, one can show after a fair bit of algebra that the advection 
and diffusion terms transform as 



D' k = %k Di + 
ah 

dh dli 



d 2 I'k n 

-A 



dhdh 



D'u, 



Di 



(9) 
(10) 



The phase-space distribution function transforms as /'(I) = Jf(l) (cf. Risken 1989) and in the new variables, the equation 
takes the standard Fokker-Planck form: 



df'(I,t) 
dt 



_d_ 



d 2 



/'(I,*)- 



(11) 



Now let I' = [E, k, cos 0). Assuming that the distribution function / is time-independent and non-zero, we may integrate 
equation (jll]) over k and cos/3. Since both of these variables have a bounded domain, the flux through their boundaries must 
vanish, leaving a single flux term: 



9(f) _ d I , d ( „ 



(12) 



where the angle brackets denote integration over k and cos /3 and sum over j denote the sum over all three variables. The 
isotropically averaged Fokker-Planck equation is then 



Of(E) 
dt 



_d_ 

dE 



d 



(D E ) is J(E) + — ((D EE ) iso f(E)) 



where {De)xso and (DEE)iso are the isotropically averaged diffusion coefficients: 



J max (-^0 

W) 



I D E 
\D EE 
where 
RE) 

and the phase-space volume is 
P(E) 



dKd(cO S 0)K^^f(E,K,0) 



Jmax(E) / dlid(cOS[3) Kf(E,K,P) 



Jmax{E) I dKd(cOS(3)K. 



(13) 



(14) 



(15) 



(16) 



Note that the standard notation in the globular cluster literature is f(E) = f(E)/P(E). 



3 NOISE MODELS 

As described in we represent distortions in the structure of halo in a biorthogonal basis. Any distortion can then be 
summarized by a set of coefficients. Because large spatial scales are most important in understanding global evolution, we can 
truncate this expansion and still recover most of the power. Internal and therefore quasi-periodic distortions contribute at a 
discrete spectrum of frequencies. Paper 1, §2 (see eqns. 21-22 in Paper 1 for the final result) derives the response of a halo to 
a point perturbation at a single frequency. Similar arguments lead to an expression for a continuous spectrum of perturbation 
frequencies. In the latter case, one computes the response of the stellar system to each frequency in the spectrum and then 
sums over all frequencies. We will begin with the development common to both cases. 

The goal is calculation of the coefficients defined by equations (§) and (|). We be gin by determining these coefficients 
for action variables and transform to (E, k, cos (5) in the end. Because orbits in the equilibrium phase space are quasi-periodic 
and representable as fixed actions and constantly advancing angles, any perturbed quantity can be represented as a Fourier 
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series in angles with coefficients depending on actions. Following Paper 1, the perturbed Hamiltonian is 
H(I,w) = ff„(I) + jri(I,w) (17) 
= H (I) + J2Hn(I)e hvr (18) 



where 



= E E E^^AOjrUCSXiiCO^W (19) 



Z=0 T7i--i j 
oo i 



= E E E F "^/ 2 ' ) r UC9X 1 S mff) / ^e 4 -* 53 [m$(u) + s jk \ 

1=0 m=-l j k 

(20) 

where 1 = h, h, h is a triple of integers, rlj(/3) and Wj^j (I) are the rotation matrices and gravitational potential transforms 
defined in Paper 1. The time dependence of the coefficients describing the response, a,j(t), is represented as its Fourier 
transform. This allows each frequency to be treated separately. The response matrix M describes the reaction of the galaxy 
to the perturbation; so the entire response is the sum of both the response and direct forcing, MI£ + 5jk- We may integrate 
the equations of motion directly to evaluate I(r + t) . Hamilton's equations yield 

dH 



h = -^- = -iY j l j Hn{l)e lv ' (21) 

OWj — ' 

i 

and therefore we have 

A/, (t + T ) = I 3 (t + r ) - I, (t ) = / + dti 3 (t ) . (22) 



The evolution of perturbed distribution function in time follows from the linearized collisionless Boltzmann equation and 
the total time derivative for a Hamiltonian system: 

i _ dfi r dfi dH dfi dHi df , , 

Analogous to the development above for Hi(t), we have 

A&w.t) =53e il w il- ^Hu(I,t) (24) 
i 

and therefore 

/ 1 (I,w,t + r)=5^e <1,w tl-^ / + dt'Hu{U) (25) 
i 

To evaluate equation (|^) we need the first and second order action moments defined by equation ^j. We assume, 
by adopting the time-asymptotic response matrix M lm in deriving that /i and Hi above, that r is larger than intrinsic 
dynamical times, consistent with the ordering of our slow and fast time scales. Previous work, including comparison to n-body 
simulations, suggests that this is a very good approximation for time scales longer than several crossing times. The response 
to the distortion induces a shift in the actions I and the overall response causes a change in the distribution function. This is 
represented in the matrix equation defined by equation ( |2c| ) for an external perturbation described by the coefficients b l ™(t). 
The fiducial stochastic variables are the coefficients fc^ m (£) themselves. 

The overall conditional probability required in equation (^) and the following development for the moments therefore has 
two contributions. First, the response of the galaxy changes the underlying distribution and consequently the probability of 
obtaining a given final state. Second, the resonant coupling changes the action of an orbit at a particular point in phase space. 
Altogether we have 



P(T, t + r\I, t) = (l + t^^J s(l'-I- £ T dt'Ai(t')\ 



(26) 



The first- and second-order moments are proportional to the square of the perturbation coefficients b 2 and are therefore 
second-order in the perturbation amplitude. With additional work, one can show that the next contributing order is propor- 
tional to b 4 and therefore relatively negligible. We will denote the ensemble average of the fluctuating coefficients which we 
will denote as angle brackets, e.g.: {b 1 ™ {tift 1 ™ fo) ■ ■ ■ b l J p% {t n )). Note that {b l ™ {ti)bf m fa)) is related to the density correlation 
function: 



{b l k m (ti)bl lm (t 2 )) = / d 3 nj d 3 r 2 YUdiAi)YU92A2)u* lm (r 1 )u[ m (r2)(p(r 1 ,t 1 )p(v 2 ,t2)}. (27) 
Assuming that the process causing fluctuations is independent of time (e.g. a stationary process) we can write 
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(p(ri,ti)p(r 3 ,t2)) = C(ri,r 2 ,ti - t 2 ). The quantity (b{ m (t 1 )b* k lm {t 2 )) = (b l k m [0)b* k lm -t 2 )) describes the correlation of 
random variables bj as a function of time. The limit ii — » t 2 gives the mean-squared fluctuation amplitude and was explored 
and compared to n-body simulations in Paper 1. 

We can now use equations (|2c|), (|2ll) and (|2^) to evaluate the moments in equation (B). After explicit substitution and 
averaging over angles, we have 



(Aljit + r)} = (^) 2 / dwi / dwa(fl 

J —oo J — oo \ 



91n/ 



91 

]T V„ 2 (tt/2, 0)rU(/?)W&£(I) Yl K m r(^) b l r m (h ) ■ 



V S j 

'"" /o |y„ 2 (7r/2,0)| 2 (i-)* / dwi / do, 2 x 



M 91 

/if rs 

t + T /*t + T 

dti y dt 2 e i(witl+W2t2) (65- m (wi)6: !m (a;2)) . (28) 

The expression for {AIj(t + r)AIu{t + r)} is nearly the same, with 1 ■ dlnf /dl in equation ( pi) ) replaced with — i^. Although 
cumbersome in appearance, all quantities in this expression are straightforwardly computed. In particular, the rotation matri- 
ces, r\ m (/3), have closed-form analytic expressions and the response matrices, M l ^(u>), have elements that can be evaluated 
by quadrature. For a specific stochastic process, all that remains is to evaluate the Fourier transform of the density correlation 
function. We will do this below for quasi-periodic and transient noise sources. 

3.1 Orbiting point mass perturbers 

For a halo of black holes, we may assume that the perturbers are point masses. In other words, the density p is a sum of delta 
functions. Expanding the distribution for a single black home in an action-angle series gives 

b l r(t) = ^y ii2 (^/2,0)r ; ' 2m (/3)W / ^ 1 J(I)e l, ' w(i) (29) 
l 

where w(t) = w + fit and after Fourier transforming in time, we find 

b i rn = Y, Y ^ n / 2 ' o y^rn(mu 2 ^(iy 1 - w °2ns(u-i-n) (so) 
i 

As in Paper 1, we assume that orbits of individual particles are uncorrelated. The particle wakes do in fact give rise to 
correlations but this is of higher order in 1/iV in the BBGKY expansion than the lowest-order effect we will consider here (cf. 
Gilbert 1969). The number density of particles at Ii,wi at time t\ and at l2,W2 at time t% is 

P(Il,Wi,ti;I 2 ,W2,t 2 ) =P(Ii,Wi)<S(Ii -I 2 )£(wi -W2 + fi(Ii)(*2-*x)) (31) 
where V(I, w) is the equilibrium particle distribution with 

N = d 3 Id 3 wV(I,w). (32) 



Direct substitution demonstrates that equation ( pi) solves the Liouville equation with the initial condition I 2 = Ii and 
w 2 = wi at t = t\. Similarly, integrating equation (pl|) over all coordinates gives N. The ensemble average (b l ™ (t\)b* 3 lm '^2)) 
is then the average of b l ™ (ti)bl lm fa) the distribution given by equation (|3l|). 

We now apply this two-particle distribution function to explicitly evaluate the ensemble average in equation The 
ensemble average here implies an average of possible distributions of point masses consistent with some underlying distribution. 
Therefore 



( Wl )&f m ' M) = 5 mm'Y ll2 (tt/2, 0)YA 2 (tt/2, 0)r\ 2m {(5)r^ {(3) x 



and 



(33) 
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^ m (ii)&: i,m (i 2 )) = {^) J <*"i J dw 2 e^ H -^ (bt m (u;i)b* s lm (u 2 )y (34) 
Integration over angles identifies 1 with 1' and therefore m = Z3 = I3 = m'. Now using equations (§|), ©, © 

to evaluate 

the average we find that the final part of expression becomes 

(i) 2 (/ ^ / dcj2ei( " ltl+ " 2t2)& ' m ^ 1 ) 6 ^ m ^ 2 )) = 
(27r) 3 ]T y d^/oWe^^"^^^^^ 



(35) 

where the integration over uj\ and W2 has become a sum over discrete frequencies denoted 1 and we have exploited the 
orthogonality of rotation matrices: 

dpsm{p)r l M r l ^((3) = ^jS u , (36) 

(Edmonds 1960). 

We now substitute this development into equation ( pj[ ) and find 

+ r)> = -^ij^ 1 ■ t£ l^frA 0)| 2 J] 2 r{ 9m 09Hm09)Wi , I ^(I)W,^ m "(l) x 

j(2^) 3 ]T f d 3 J/ o (I)|y„ 2 (7r/2,0)]VU(/^^ * 



1 

t + T ^"t + T 

ii.n(i)(ti-t 2 ) 



dh / d£ 2 e"'" wltl - t2; ^ . (37) 



As described generally above, the term in {■} is the temporal correlation function for the response to the point masses 
fluctuations and only depends on the time difference t\ — £2. The correlation is finite and lasts some order-unity number 
of dynamical times. Our diffusion calculation is in the regime r S> 1/fi. We may change the double time integration from 
variables £1, £2 to T = (£1 +£2)/2, r = £1 — £2. The integral over r gives a delta function (5(1-0(1)) and the integral over T gives 
t. The delta function implies that only orbits with commensurate frequencies will give rise to secular changes. Physically, the 
disturbance must present an asymmetric force distribution to cause a secular change in the actions of dark matter orbits. If 
the orbital frequencies are not commensurate, the long-term average of the perturbing force will be axisymmetric. For orbiting 
black holes, the ratio ranges from 1 to 2 as the halo potential varies from that for point mass to homogeneous core. 

Therefore for most systems, no commensurabilities are available for harmonic orders I = 1, 2. In other words, a resonant I = 1 
orbit will look like a Keplerian orbit and a resonant I = 2 orbit will look like a bisymmetric oval (e.g. stationary bar orbit) 
and neither exist in any significant measure in most extended stellar systems. 

Returning to equation (Q), we can now read off our diffusion coefficients in action variables: 

Df\l,t) = lim ££g±I» 

= ' w |y " a(V2 ' 0)|2 £ E ••i™(ffl»-t(fflw i t(i)< i ;(i) x 

fj.b' rs 

|(2^) 3 ^] j d 3 If o (I)\Y ll2 (7T/2,0)\ 2 r l hm (P)r l l2m (f3)x 



(38) 



D (,) (M) = lim (A7 j( £ + r)A7 fc (£ + r)) 
J t— >o 2r 
Ijlk 



[if rs 

|(2^) 3 ^] J d 3 If o {I)\Yu 2 {7v/2,0)\ 2 r l l2m {P)r l l2m (P)x 

<:(i)<;(p^(i ■ m)M: i rd ■ m)^s (i . n(i))} . 



(39) 
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Note that the limit r — > is taken in the sense that r is small compared to the evolutionary time scale due to the fluctuations 
but large compared to the dynamical time. The time dependence in the diffusion coefficients reminds us that the underlying 
equilibrium distribution /o(I) changes on an evolutionary time scale but, for the purposes of computation, is held fixed on a 
dynamical time scale. The integrals may be simplified by noting that d 3 I — dEdJ Jd(cos /3)/fii (E, J), we can do the integral 
in P using the orthogonality of the rotation matrices as previously described. For a given background distribution function 
/ (I), the term in curly brackets can be computed once and for all since they are independent of the local value of the actions. 



3.2 Transient processes 

In the previous application for orbiting point masses, we saw that only harmonics with I > 3 were effective at driving evolution. 
Here, we describe the results of bombarding the galaxy isotropically with bits of mass. This is an idealization of interactions 
during the epoch of galaxy evolution. The perturbations are shots, not quasi-periodic, and therefore lead to a continuous 
spectrum of perturbation frequencies. This case and the previously considered point-mass case represent two extremes. For 
example, a decaying orbit will have a set of broad peaks and a low-frequency continuum. 

The expansion coefficients for the perturbation are straightforward for the mass inside of the galaxy: 

bT(t) = J ' d\YUe,<t>)u* lm {r)5{r -r(t)) = YUO(t),rtt)K lrn (r(t)). (40) 

If the mass is outside the galaxy, we use the multipole expansion with the density rather than the potential member of the 
biorthogonal pair to evaluate the coefficients: 

b\ m {t) = ~ ( d 3 rYC m (9,m lm (r) ' 



4tt J VJ \r-r(t)\ 

« i — i > 



4tt 

1 



Yr m (e(t),<i>(t))YC m (e,<i>) 



2/ + f l£,(0(t),*(t)) J drr 2 d* 1 ™^)^-^ (41) 

The Fourier transform needed for equation (|i^) is most easily done assuming the perturber is in the equatorial plane. 
We denote the Fourier transform of b 1 ™ 1 ^) as b^ (a;). In practice, we will perform the transform numerically by FFT. Note 
that equations ( fib] ) and ( |4l| ) describe the time dependence in coefficients for any trajectory r(t),8(t),<j)(t). A cloud of points 
and more generally any phase-space distortion yielding bl m (t) or coefficients at discrete times from an n-body simulation 
are possible input to the FFT. Now, we can evaluate the final line in equation ( |28| ) by changing coordinates from ti,ts to 
T = (ti +t 2 )/2,r =ti- t 2 . We have 

i(a>lti+w 2 t 3 ) _ e »("i[T+r/2]+o ;2 [T-T/2]) _ e i(u J1 +uj 2 )T e i(u J1 ~u J2 )T/2 

Using this in the last line of equation ( jj^ ) we find 

/t+r rt + r/2 

di 2 e l( " ltl+ " 2t2) (^ m (^i)b:' m M) = J dTAnSfa - ^ 2 )e l2 " T (b^O^M) 

= 4tt5(uji — U2)e^ T ——— (V™ (tvi)b* lm (lu2)) ■ 

(42) 

where u> = lj 1 — ll>2- In deriving the second equality, we note that the bombardment must occur between t and t + r and use 
the shift properties of the Fourier transform to write (b l ™ (wi)b* im (w 2 )) = e~ 2bJt (b l r m (u!i)b* s lm (oj 2 )) where the transform b{to) 
denotes the transform of an event centred about the temporal origin. In the limit r — > 0, this expression becomes 4nS(cui — u>2)t. 
Substituting, this back into equation (E8|) we can perform one of the uj integrals straight away. After rearranging we have 



(AIj(t + r)> = -J^IA ■ ^ \Yu 2 (tt/2, 0)| 2 ^ £ rUM^M^^lXJfl) x 



liv r 
i ? >J5 



4tt / «^HA<*H(f( Wl )i; im (^)). (43) 



The expression for {AIj(t + r)AIk(t + r)) follows by analogy with the equations ( |2q ) and (g9|) in §3.1. 

Symmetry suggests the choice of perturbing orbits on the equatorial plane should not effect the final results. This can 
be explicitly demonstrated explicitly using the rotational properties of the spherical harmonics. Let 1Z l mrn i (a, /?, 7) be the 
rotation matrix with the Euler angles a, /3, 7. Then, because 



1 
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where the primed coordinates refer to the rotated coordinate system, we have 
i 

m' '= — I 

For an isotropic spherical system M 1 ™{uj) is independent of m. We can exploit this and the rotational properties of b l ™(ui) 
to simplify the computation of (b l ™ (ijj\)b* lm (a^)) . We may express bj. m ((j;) in any convenient coordinate system and use the 
rotation matrices to rotate this to any orientation: 

(b l r m (^)b: lm (u; 2 )) = / C(",ft7)t'K) Yl Rl m >>(a,P,l)K lm ''(Lo 2 )\ (45) 

\ m ; — — I m ff =—l I 

Summing over all values m for a given I we have 

/ b'T(cu 1 )b: lm (uj 2 )\ = {Y1 &^V)«' m ' (<■*)) (46) 

\m=-i / \m' = -l I 

having used the orthogonality of rotation matrices. Finally, we assume that the ensemble average includes random events 
from all directions and therefore only the same event will be correlated in the computation of ^fcj" 1 {u)\)b* s lrn '(u^)) • 



4 SUMMARY 

This paper presents a general equation for noise-driven evolution which incorporates the full self-gravitating response of a 
stochastic process in the perturbation limit. The motivating question is same as that of violent relaxation: what is the long- 
term response of a stellar system to a fluctuating potential? By working in the perturbation limit, the linearity guarantees that 
the process is Markovian and therefore the response of the stellar system to repeated stochastic events is naturally treated 
using the matrix-method and the statistic methods for handling general stochastic differential equations. 

In general, expansions of integral equations for stochastic processes yield an infinite number of terms. A general theorem 
demonstrates that the truncation at quadratic order is consistent for a Markov process and the resulting evolution equation 
takes the Fokker-Planck form0. Evaluation of the Fokker-Planck coefficients requires the specifying the response of stellar 
system individual events in the stochastic process. 

We explicitly develop techniques for two general situations: quasi-periodic perturbers and transient perturbers. The former 
case is motivated by halo of super-massive black holes (e.g. Lacey & Ostriker 1985). The latter case includes almost everything 
else; for example, unbound dwarf encounters (fly-by), orbiting substructure decaying due to dynamical friction, disrupting 
dwarfs, and mixing tidal debris. There is no practical constraint on deriving the Fokker-Planck coefficients for transient noise 
but that the process can be represented as an expansion in some biorthogonal basis. For dwarfs on decaying or unbound 
orbits, this is particularly easy and can be done by quadrature. Alternatively, one may construct an n-body simulation using 
the expansion method and let the simulation produce the time series of coefficients directly. 

The companion paper (Paper 2) applies the apparatus described here to investigate the evolution of haloes during the 
noisy epoch of galaxy formation. We find that both unbound encounters and decaying substructure drives the halo profile to 
a self-similar form similar to those found recently in cosmological simulations. There are a number of interesting applications. 
For example, the distribution of binary semi-major axes a in the field star populations is proportional to a -1 over several orders 
of magnitude. Preliminary work suggests that this may be explained as the result of fluctuations from the noisy environment 
found in proto-stellar clusters and molecular clouds. Other possible applications include effect of transient bar formation and 
spiral structure on overall galactic structure in the 10 billion years since formation. 
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